Usenet 1993 July/InfoMagic USENET CD-ROM July 1993.ISO/sources/misc/volume13/gmcalc/part08

Text File | 1990-06-05 | 57.1 KB | 1,896 lines

Newsgroups: comp.sources.misc
From: [e-mail] (David Gillespie)
Subject: v13i034: Emacs Calculator 1.01, part 08/19
Sender: [e-mail] (Brandon S. Allbery - comp.sources.misc)

Posting-number: Volume 13, Issue 34
Submitted-by: [e-mail] (David Gillespie)
Archive-name: gmcalc/part08

---- Cut Here and unpack ----
#!/bin/sh
# this is part 8 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc-ext.el continued
#
CurArch=8
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
     exit 1; fi
( read Scheck
  if test "$Scheck" != $CurArch
  then echo "Please unpack part $Scheck next!"
       exit 1;
  else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc-ext.el"
sed 's/^X//' << 'SHAR_EOF' >> calc-ext.el
X                   (- (+ (math-numdigs (nth 1 a))
X                     (nth 2 a))
X                      (1+ tol)))))
X    ((not (eq (car tol) 'float))
X     (if (Math-realp tol)
X         (math-to-fraction a (math-float tol))
X       (math-reject-arg tol 'realp)))
X    ((Math-negp tol)
X     (math-to-fraction a (math-neg tol)))
X    ((Math-zerop tol)
X     (math-to-fraction a 0))
X    ((not (math-lessp-float tol '(float 1 0)))
X     (math-trunc a))
X    ((Math-zerop a)
X     0)
X    (t
X     (let ((cfrac (math-continued-fraction a tol))
X           (calc-prefer-frac t))
X       (math-eval-continued-fraction cfrac))))
X)
X(fset 'calcFunc-frac (symbol-function 'math-to-fraction))
X
X(defun math-continued-fraction (a tol)
X  (let ((calc-internal-prec (+ calc-internal-prec 2)))
X    (let ((cfrac nil)
X      (aa a)
X      (calc-prefer-frac nil)
X      int)
X      (while (or (null cfrac)
X         (and (not (Math-zerop aa))
X              (not (math-lessp-float
X                (math-abs
X                 (math-sub a
X                       (let ((f (math-eval-continued-fraction
X                         cfrac)))
X                     (math-working "Fractionalize" f)
X                     f)))
X                tol))))
X    (setq int (math-trunc aa)
X          aa (math-sub aa int)
X          cfrac (cons int cfrac))
X    (or (Math-zerop aa)
X        (setq aa (math-div 1 aa))))
X      cfrac))
X)
X
X(defun math-eval-continued-fraction (cf)
X  (let ((n (car cf))
X    (d 1)
X    temp)
X    (while (setq cf (cdr cf))
X      (setq temp (math-add (math-mul (car cf) n) d)
X        d n
X        n temp))
X    (math-div n d))
X)
X
X
X(defun math-clean-number (a &optional prec)   ; [X X S] [Public]
X  (if prec
X      (cond ((Math-messy-integerp prec)
X         (math-clean-number a (math-trunc prec)))
X        ((or (not (integerp prec))
X         (< prec 3))
X         (calc-record-why "Precision must be an integer 3 or above")
X         (list 'calcFunc-clean a prec))
X        ((not (Math-objvecp a))
X         (list 'calcFunc-clean a prec))
X        (t (let ((calc-internal-prec prec))
X         (math-clean-number (math-normalize a)))))
X    (cond ((eq (car-safe a) 'polar)
X       (let ((theta (math-mod (nth 2 a)
X                  (if (eq calc-angle-mode 'rad)
X                      (math-two-pi)
X                    360))))
X         (math-neg
X          (math-neg
X           (math-normalize
X        (list 'polar (nth 1 a) theta))))))
X      ((Math-vectorp a)
X       (math-map-vec 'math-clean-number a))
X      ((Math-objectp a) a)
X      (t (list 'calcFunc-clean a))))
X)
X(fset 'calcFunc-clean (symbol-function 'math-clean-number))
X
X
X
X
X;;;; Logical operations.
X
X(defun calcFunc-eq (a b)
X  (let ((res (math-compare a b)))
X    (if (= res 0)
X    1
X      (if (= res 2)
X      (list 'calcFunc-eq a b)
X    0)))
X)
X
X(defun calcFunc-neq (a b)
X  (let ((res (math-compare a b)))
X    (if (= res 0)
X    0
X      (if (= res 2)
X      (list 'calcFunc-neq a b)
X    1)))
X)
X
X(defun calcFunc-lt (a b)
X  (let ((res (math-compare a b)))
X    (if (= res -1)
X    1
X      (if (= res 2)
X      (list 'calcFunc-lt a b)
X    0)))
X)
X
X(defun calcFunc-gt (a b)
X  (let ((res (math-compare a b)))
X    (if (= res 1)
X    1
X      (if (= res 2)
X      (list 'calcFunc-gt a b)
X    0)))
X)
X
X(defun calcFunc-leq (a b)
X  (let ((res (math-compare a b)))
X    (if (= res 1)
X    0
X      (if (= res 2)
X      (list 'calcFunc-leq a b)
X    1)))
X)
X
X(defun calcFunc-geq (a b)
X  (let ((res (math-compare a b)))
X    (if (= res -1)
X    0
X      (if (= res 2)
X      (list 'calcFunc-geq a b)
X    1)))
X)
X
X(defun calcFunc-land (a b)
X  (cond ((Math-zerop a)
X     a)
X    ((Math-zerop b)
X     b)
X    ((Math-numberp a)
X     b)
X    ((Math-numberp b)
X     a)
X    (t (list 'calcFunc-land a b)))
X)
X
X(defun calcFunc-lor (a b)
X  (cond ((Math-zerop a)
X     b)
X    ((Math-zerop b)
X     a)
X    ((Math-numberp a)
X     a)
X    ((Math-numberp b)
X     b)
X    (t (list 'calcFunc-lor a b)))
X)
X
X(defun calcFunc-lnot (a)
X  (if (Math-zerop a)
X      1
X    (if (Math-numberp a)
X    0
X      (list 'calcFunc-lnot a)))
X)
X
X(defun calcFunc-if (c e1 e2)
X  (if (Math-zerop c)
X      e2
X    (if (Math-numberp c)
X    e1
X      (list 'calcFunc-if c e1 e2)))
X)
X
X(defun math-normalize-logical-op (a)
X  (or (and (eq (car a) 'calcFunc-if)
X       (= (length a) 4)
X       (let ((a1 (math-normalize (nth 1 a))))
X         (if (Math-zerop a1)
X         (math-normalize (nth 3 a))
X           (if (Math-numberp a1)
X           (math-normalize (nth 2 a))
X         (list 'calcFunc-if a1 (nth 2 a) (nth 3 a))))))
X      a)
X)
X
X(defun calcFunc-in (a b)
X  (or (and (eq (car-safe b) 'vec)
X       (let ((bb b))
X         (while (and (setq bb (cdr bb))
X             (not (if (memq (car-safe (car bb)) '(vec intv))
X                  (eq (calcFunc-in a (car bb)) 1)
X                (math-equal a (car bb))))))
X         (if bb 1 (and (math-constp a) (math-constp bb) 0))))
X      (and (eq (car-safe b) 'intv)
X       (let ((res (math-compare a (nth 2 b))))
X         (cond ((= res -1)
X            0)
X           ((= res 0)
X            (if (memq (nth 1 b) '(2 3)) 1 0))
X           ((/= res 1)
X            nil)
X           ((= (setq res (math-compare a (nth 3 b))) 1)
X            0)
X           ((= res 0)
X            (if (memq (nth 1 b) '(1 3)) 1 0))
X           ((/= res -1)
X            nil)
X           (t 1))))
X      (and (math-equal a b)
X       1)
X      (and (math-constp a) (math-constp b)
X       0)
X      (list 'calcFunc-in a b))
X)
X
X(defun calcFunc-typeof (a)
X  (cond ((Math-integerp a) 1)
X    ((eq (car a) 'frac) 2)
X    ((eq (car a) 'float) 3)
X    ((eq (car a) 'hms) 4)
X    ((eq (car a) 'cplx) 5)
X    ((eq (car a) 'polar) 6)
X    ((eq (car a) 'sdev) 7)
X    ((eq (car a) 'intv) 8)
X    ((eq (car a) 'mod) 9)
X    ((eq (car a) 'var) 100)
X    ((eq (car a) 'vec) (if (math-matrixp a) 102 101))
X    (t (let ((func (assq (car a) '( ( + . calcFunc-add )
X                    ( - . calcFunc-sub )
X                    ( * . calcFunc-mul )
X                    ( / . calcFunc-div )
X                    ( ^ . calcFunc-pow )
X                    ( % . calcFunc-mod )
X                    ( neg . calcFunc-neg )
X                    ( | . calcFunc-vconcat ) ))))
X         (setq func (if func (cdr func) (car a)))
X         (math-calcFunc-to-var func))))
X)
X
X(defun calcFunc-integer (a)
X  (if (Math-integerp a)
X      1
X    (if (Math-objvecp a)
X    0
X      (list 'calcFunc-integer a)))
X)
X
X(defun calcFunc-real (a)
X  (if (Math-realp a)
X      1
X    (if (Math-objvecp a)
X    0
X      (list 'calcFunc-real a)))
X)
X
X(defun calcFunc-constant (a)
X  (if (math-constp a)
X      1
X    (if (Math-objvecp a)
X    0
X      (list 'calcFunc-constant a)))
X)
X
X(defun calcFunc-refers (a b)
X  (if (math-expr-contains a b)
X      1
X    (if (eq (car-safe a) 'var)
X    (list 'calcFunc-refers a b)
X      0))
X)
X
X
X
X
X;;;; Complex numbers.
X
X(defun math-to-polar (a)   ; [C N] [Public]
X  (cond ((Math-vectorp a)
X     (math-map-vec 'math-to-polar a))
X    ((Math-realp a) a)
X    ((Math-numberp a)
X     (math-normalize (math-polar a)))
X    (t (list 'calcFunc-polar)))
X)
X(fset 'calcFunc-polar (symbol-function 'math-to-polar))
X
X(defun math-to-rectangular (a)   ; [N N] [Public]
X  (cond ((Math-vectorp a)
X     (math-map-vec 'math-to-rectangular a))
X    ((Math-realp a) a)
X    ((Math-numberp a)
X     (math-normalize (math-complex a)))
X    (t (list 'calcFunc-rect)))
X)
X(fset 'calcFunc-rect (symbol-function 'math-to-rectangular))
X
X;;; Compute the complex conjugate of A.  [O O] [Public]
X(defun math-conj (a)
X  (cond ((math-real-objectp a)
X     a)
X    ((eq (car a) 'cplx)
X     (list 'cplx (nth 1 a) (math-neg (nth 2 a))))
X    ((eq (car a) 'polar)
X     (list 'polar (nth 1 a) (math-neg (nth 2 a))))
X    ((eq (car a) 'vec)
X     (math-map-vec 'math-conj a))
X    ((eq (car a) 'calcFunc-conj)
X     (nth 1 a))
X    (t (calc-record-why 'numberp a)
X       (list 'calcFunc-conj a)))
X)
X(fset 'calcFunc-conj (symbol-function 'math-conj))
X
X;;; Compute the absolute value squared of A.  [F N] [Public]
X(defun math-abssqr (a)
X  (cond ((Math-realp a)
X     (math-sqr a))
X    ((eq (car a) 'cplx)
X     (math-add (math-sqr (nth 1 a)) (math-sqr (nth 2 a))))
X    ((eq (car a) 'polar)
X     (math-sqr (nth 1 a)))
X    ((and (memq (car a) '(sdev intv)) (math-constp a))
X     (math-sqr (math-abs a)))
X    ((eq (car a) 'vec)
X     (math-reduce-vec 'math-add (math-map-vec 'math-abssqr a)))
X    (t (calc-record-why 'numvecp a)
X       (list 'calcFunc-abssqr a)))
X)
X(fset 'calcFunc-abssqr (symbol-function 'math-abssqr))
X
X;;; Compute the complex argument of A.  [F N] [Public]
X(defun math-cplx-arg (a)
X  (cond ((Math-anglep a)
X     (if (math-negp a) (math-half-circle nil) 0))
X    ((eq (car-safe a) 'cplx)
X     (math-arctan2 (nth 2 a) (nth 1 a)))
X    ((eq (car-safe a) 'polar)
X     (nth 2 a))
X    ((eq (car a) 'vec)
X     (math-map-vec 'math-cplx-arg a))
X    (t (calc-record-why 'numvecp a)
X       (list 'calcFunc-arg a)))
X)
X(fset 'calcFunc-arg (symbol-function 'math-cplx-arg))
X
X;;; Extract the real or complex part of a complex number.  [R N] [Public]
X;;; Also extracts the real part of a modulo form.
X(defun math-real-part (a)
X  (cond ((memq (car-safe a) '(mod sdev))
X     (nth 1 a))
X    ((math-real-objectp a) a)
X    ((eq (car a) 'cplx)
X     (nth 1 a))
X    ((eq (car a) 'polar)
X     (math-mul (nth 1 a) (math-cos (nth 2 a))))
X    ((eq (car a) 'vec)
X     (math-map-vec 'math-real-part a))
X    (t (calc-record-why 'numberp a)
X       (list 'calcFunc-re a)))
X)
X(fset 'calcFunc-re (symbol-function 'math-real-part))
X
X(defun math-imag-part (a)
X  (cond ((math-real-objectp a)
X     (if (math-floatp a) '(float 0 0) 0))
X    ((eq (car a) 'cplx)
X     (nth 2 a))
X    ((eq (car a) 'polar)
X     (math-mul (nth 1 a) (math-sin (nth 2 a))))
X    ((eq (car a) 'vec)
X     (math-map-vec 'math-imag-part a))
X    (t (calc-record-why 'numberp a)
X       (list 'calcFunc-im a)))
X)
X(fset 'calcFunc-im (symbol-function 'math-imag-part))
X
X
X
X;;;; Transcendental functions.
X
X;;; All of these functions are defined on the complex plane.
X;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
X
X;;; Most functions increase calc-internal-prec by 2 digits, then round
X;;; down afterward.  "-raw" functions use the current precision, require
X;;; their arguments to be in float (or complex float) format, and always
X;;; work in radians (where applicable).
X
X(defun math-to-radians (a)   ; [N N]
X  (cond ((eq (car-safe a) 'hms)
X     (math-from-hms a 'rad))
X    ((memq calc-angle-mode '(deg hms))
X     (math-mul a (math-pi-over-180)))
X    (t a))
X)
X
X(defun math-from-radians (a)   ; [N N]
X  (cond ((eq calc-angle-mode 'deg)
X     (if (math-constp a)
X         (math-div a (math-pi-over-180))
X       (list 'calcFunc-deg a)))
X    ((eq calc-angle-mode 'hms)
X     (math-to-hms a 'rad))
X    (t a))
X)
X
X(defun math-to-radians-2 (a)   ; [N N]
X  (cond ((eq (car-safe a) 'hms)
X     (math-from-hms a 'rad))
X    ((memq calc-angle-mode '(deg hms))
X     (if calc-symbolic-mode
X         (math-div (math-mul a '(var pi var-pi)) 180)
X       (math-mul a (math-pi-over-180))))
X    (t a))
X)
X
X(defun math-from-radians-2 (a)   ; [N N]
X  (cond ((memq calc-angle-mode '(deg hms))
X     (if calc-symbolic-mode
X         (math-div (math-mul 180 a) '(var pi var-pi))
X       (math-div a (math-pi-over-180))))
X    (t a))
X)
X
X
X
X;;; Sine, cosine, and tangent.
X
X(defun math-sin (x)   ; [N N] [Public]
X  (cond ((Math-scalarp x)
X     (math-with-extra-prec 2
X       (math-sin-raw (math-to-radians (math-float x)))))
X    ((eq (car x) 'sdev)
X     (if (math-constp x)
X         (math-with-extra-prec 2
X           (let* ((xx (math-to-radians (math-float (nth 1 x))))
X              (xs (math-to-radians (math-float (nth 2 x))))
X              (sc (math-sin-cos-raw xx)))
X         (math-make-sdev (car sc)
X                 (math-mul xs (math-abs (cdr sc))))))
X       (math-make-sdev (math-sin (nth 1 x))
X               (math-mul (nth 2 x)
X                     (math-abs (math-cos (nth 1 x)))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-cos (math-sub x (math-quarter-circle nil))))
X    (t (calc-record-why 'scalarp x)
X       (list 'calcFunc-sin x)))
X)
X(fset 'calcFunc-sin (symbol-function 'math-sin))
X
X(defun math-cos (x)   ; [N N] [Public]
X  (cond ((Math-scalarp x)
X     (math-with-extra-prec 2
X       (math-cos-raw (math-to-radians (math-float x)))))
X    ((eq (car x) 'sdev)
X     (if (math-constp x)
X         (math-with-extra-prec 2
X           (let* ((xx (math-to-radians (math-float (nth 1 x))))
X              (xs (math-to-radians (math-float (nth 2 x))))
X              (sc (math-sin-cos-raw xx)))
X         (math-make-sdev (cdr sc)
X                 (math-mul xs (math-abs (car sc))))))
X       (math-make-sdev (math-cos (nth 1 x))
X               (math-mul (nth 2 x)
X                     (math-abs (math-sin (nth 1 x)))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-with-extra-prec 2
X       (let* ((xx (math-to-radians (math-float x)))
X          (na (math-floor (math-div (nth 2 xx) (math-pi))))
X          (nb (math-floor (math-div (nth 3 xx) (math-pi))))
X          (span (math-sub nb na)))
X         (if (memq span '(0 1))
X         (let ((int (math-sort-intv (nth 1 x)
X                        (math-cos-raw (nth 2 xx))
X                        (math-cos-raw (nth 3 xx)))))
X           (if (eq span 1)
X               (if (math-evenp na)
X               (math-make-intv (logior (nth 1 x) 2)
X                       -1
X                       (nth 3 int))
X             (math-make-intv (logior (nth 1 x) 1)
X                     (nth 2 int)
X                     1))
X             int))
X           (list 'intv 3 -1 1)))))
X    (t (calc-record-why 'scalarp x)
X       (list 'calcFunc-cos x)))
X)
X(fset 'calcFunc-cos (symbol-function 'math-cos))
X
X(defun math-sin-cos (x)   ; [V N] [Public]
X  (if (Math-scalarp x)
X      (math-with-extra-prec 2
X    (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
X      (list 'vec (cdr sc) (car sc))))    ; the vector [cos, sin]
X    (list 'vec (math-sin x) (math-cos x)))
X)
X(fset 'calcFunc-sincos (symbol-function 'math-sin-cos))
X
X(defun math-tan (x)   ; [N N] [Public]
X  (cond ((Math-scalarp x)
X     (math-with-extra-prec 2
X       (math-tan-raw (math-to-radians (math-float x)))))
X    ((eq (car x) 'sdev)
X     (if (math-constp x)
X         (math-with-extra-prec 2
X           (let* ((xx (math-to-radians (math-float (nth 1 x))))
X              (xs (math-to-radians (math-float (nth 2 x))))
X              (sc (math-sin-cos-raw xx)))
X         (if (math-zerop (cdr sc))
X             (progn
X               (calc-record-why "Division by zero")
X               (list 'calcFunc-tan x))
X           (math-make-sdev (math-div-float (car sc) (cdr sc))
X                   (math-div-float xs (math-sqr (cdr sc)))))))
X       (math-make-sdev (math-tan (nth 1 x))
X               (math-div (nth 2 x)
X                     (math-sqr (math-cos (nth 1 x)))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (or (math-with-extra-prec 2
X           (let* ((xx (math-to-radians (math-float x)))
X              (na (math-floor (math-div (math-sub (nth 2 xx)
X                              (math-pi-over-2))
X                        (math-pi))))
X              (nb (math-floor (math-div (math-sub (nth 3 xx)
X                              (math-pi-over-2))
X                        (math-pi)))))
X         (and (equal na nb)
X              (math-sort-intv (nth 1 x)
X                      (math-tan-raw (nth 2 xx))
X                      (math-tan-raw (nth 3 xx))))))
X         (progn
X           (calc-record-why "Infinite interval" x)
X           (list 'calcFunc-tan x))))
X    (t (calc-record-why 'scalarp x)
X       (list 'calcFunc-tan x)))
X)
X(fset 'calcFunc-tan (symbol-function 'math-tan))
X
X(defun math-sin-raw (x)   ; [N N]
X  (cond ((eq (car-safe x) 'cplx)
X     (let* ((expx (math-exp-raw (nth 2 x)))
X        (expmx (math-div-float '(float 1 0) expx))
X        (sc (math-sin-cos-raw (nth 1 x))))
X       (list 'cplx
X         (math-mul-float (car sc)
X                 (math-mul-float (math-sub expx expmx)
X                         '(float 5 -1)))
X         (math-mul-float (cdr sc)
X                 (math-mul-float (math-add-float expx expmx)
X                         '(float 5 -1))))))
X    ((eq (car-safe x) 'polar)
X     (math-polar (math-sin-raw (math-complex x))))
X    ((Math-integer-negp (nth 1 x))
X     (math-neg-float (math-sin-raw (math-neg-float x))))
X    ((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
X     (math-sin-raw (math-mod x (math-two-pi))))
X    (t (math-sin-raw-2 x x)))
X)
X
X(defun math-cos-raw (x)   ; [N N]
X  (if (eq (car-safe x) 'polar)
X      (math-polar (math-cos-raw (math-complex x)))
X    (math-sin-raw (math-sub-float (math-pi-over-2) x)))
X)
X
X;;; This could use a smarter method:  Reduce x as in math-sin-raw, then
X;;;   compute either sin(x) or cos(x), whichever is smaller, and compute
X;;;   the other using the identity sin(x)^2 + cos(x)^2 = 1.
X(defun math-sin-cos-raw (x)   ; [F.F F]  (result is (sin x . cos x))
X  (cons (math-sin-raw x) (math-cos-raw x))
X)
X
X(defun math-tan-raw (x)   ; [N N]
X  (cond ((eq (car-safe x) 'cplx)
X     (let* ((x (math-mul-float x '(float 2 0)))
X        (expx (math-exp-raw (nth 2 x)))
X        (expmx (math-div-float '(float 1 0) expx))
X        (sc (math-sin-cos-raw (nth 1 x)))
X        (d (math-add-float (cdr sc)
X                   (math-mul-float (math-add-float expx expmx)
X                           '(float 5 -1)))))
X       (and (not (eq (nth 1 d) 0))
X        (list 'cplx
X              (math-div-float (car sc) d)
X              (math-div-float (math-mul-float (math-add expx expmx)
X                              '(float 5 -1)) d)))))
X    ((eq (car-safe x) 'polar)
X     (math-polar (math-tan-raw (math-complex x))))
X    (t
X     (let ((sc (math-sin-cos-raw x)))
X       (if (eq (nth 1 (cdr sc)) 0)
X           (math-reject-arg x "Division by zero")
X         (math-div-float (car sc) (cdr sc))))))
X)
X
X(defun math-sin-raw-2 (x orgx)   ; This avoids poss of inf recursion.  [F F]
X  (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
X    (cond ((Math-integer-negp (nth 1 xmpo2))
X       (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
X                       orgx)))
X      ((math-lessp-float (math-pi-over-4) x)
X       (math-cos-raw-2 xmpo2 orgx))
X      ((math-lessp-float x (math-neg (math-pi-over-4)))
X       (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
X      ((math-nearly-zerop-float x orgx) '(float 0 0))
X      (calc-symbolic-mode (signal 'inexact-result nil))
X      (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))
X)
X
X(defun math-cos-raw-2 (x orgx)   ; [F F]
X  (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
X    (calc-symbolic-mode (signal 'inexact-result nil))
X    (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
X         (math-sin-series
X          (math-add-float '(float 1 0)
X                  (math-mul-float xnegsqr '(float 5 -1)))
X          24 5 xnegsqr xnegsqr))))
X)
X
X(defun math-sin-series (sum nfac n x xnegsqr)
X  (math-working "sin" sum)
X  (let* ((nextx (math-mul-float x xnegsqr))
X     (nextsum (math-add-float sum (math-div-float nextx
X                              (math-float nfac)))))
X    (if (math-nearly-equal-float sum nextsum)
X    sum
X      (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
X               (+ n 2) nextx xnegsqr)))
X)
X
X
X;;; Inverse sine, cosine, tangent.
X
X(defun math-arcsin (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (math-from-radians (math-arcsin-raw (math-float x)))))
X    ((and (eq (car x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (not (Math-lessp 1 (math-abs (nth 1 x))))))
X     (math-make-sdev (math-arcsin (nth 1 x))
X             (math-from-radians
X              (math-div (nth 2 x)
X                    (math-sqrt
X                     (math-sub 1 (math-sqr (nth 1 x))))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-arcsin (nth 2 x))
X             (math-arcsin (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-arcsin x)))
X)
X(fset 'calcFunc-arcsin (symbol-function 'math-arcsin))
X
X(defun math-arccos (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (math-from-radians (math-arccos-raw (math-float x)))))
X    ((and (eq (car x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (not (Math-lessp 1 (math-abs (nth 1 x))))))
X     (math-make-sdev (math-arccos (nth 1 x))
X             (math-from-radians
X              (math-div (nth 2 x)
X                    (math-sqrt
X                     (math-sub 1 (math-sqr (nth 1 x))))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-arccos (nth 2 x))
X             (math-arccos (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-arccos x)))
X)
X(fset 'calcFunc-arccos (symbol-function 'math-arccos))
X
X(defun math-arctan (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (math-from-radians (math-arctan-raw (math-float x)))))
X    ((eq (car x) 'sdev)
X     (math-make-sdev (math-arctan (nth 1 x))
X             (math-from-radians
X              (math-div (nth 2 x)
X                    (math-add 1 (math-sqr (nth 1 x)))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-arctan (nth 2 x))
X             (math-arctan (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-arctan x)))
X)
X(fset 'calcFunc-arctan (symbol-function 'math-arctan))
X
X(defun math-arcsin-raw (x)   ; [N N]
X  (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
X    (if (or (memq (car-safe x) '(cplx polar))
X        (memq (car-safe a) '(cplx polar)))
X    (math-with-extra-prec 2   ; use extra precision for difficult case
X      (math-mul '(cplx 0 -1)
X            (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
X      (math-arctan2-raw x a)))
X)
X
X(defun math-arccos-raw (x)   ; [N N]
X  (math-sub (math-pi-over-2) (math-arcsin-raw x))
X)
X
X(defun math-arctan-raw (x)   ; [N N]
X  (cond ((memq (car-safe x) '(cplx polar))
X     (math-with-extra-prec 2   ; extra-extra
X       (math-mul '(cplx 0 -1)
X             (math-ln-raw (math-mul
X                   (math-add 1 (math-mul '(cplx 0 1) x))
X                   (math-sqrt-raw
X                    (math-div 1 (math-add
X                         1 (math-sqr x)))))))))
X    ((Math-integer-negp (nth 1 x))
X     (math-neg-float (math-arctan-raw (math-neg-float x))))
X    ((math-zerop x) x)
X    ((math-equal-int x 1) (math-pi-over-4))
X    ((math-equal-int x -1) (math-neg (math-pi-over-4)))
X    (calc-symbolic-mode (signal 'inexact-result nil))
X    ((math-lessp-float '(float 414214 -6) x)  ; if x > sqrt(2) - 1, reduce
X     (if (math-lessp-float '(float 1 0) x)
X         (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
X                 (math-arctan-raw (math-div-float '(float 1 0) x)))
X       (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
X               (math-arctan-raw (math-div-float
X                         (math-sub-float '(float 1 0) x)
X                         (math-add-float '(float 1 0)
X                                 x))))))
X    (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))
X)
X
X(defun math-arctan-series (sum n x xnegsqr)
X  (math-working "arctan" sum)
X  (let* ((nextx (math-mul-float x xnegsqr))
X     (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X    (if (math-nearly-equal-float sum nextsum)
X    sum
X      (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))
X)
X
X(defun math-arctan2 (y x)   ; [F R R] [Public]
X  (if (Math-anglep y)
X      (if (Math-anglep x)
X      (math-with-extra-prec 2
X        (math-from-radians (math-arctan2-raw (math-float y)
X                         (math-float x))))
X    (calc-record-why 'anglep x)
X    (list 'calcFunc-arctan2 y x))
X    (calc-record-why 'anglep y)
X    (list 'calcFunc-arctan2 y x))
X)
X(fset 'calcFunc-arctan2 (symbol-function 'math-arctan2))
X
X(defun math-arctan2-raw (y x)   ; [F R R]
X  (cond ((math-zerop y)
X     (if (math-negp x) (math-pi) 0))
X    ((math-zerop x)
X     (if (math-posp y)
X         (math-pi-over-2)
X       (math-neg (math-pi-over-2))))
X    ((math-posp x)
X     (math-arctan-raw (math-div-float y x)))
X    ((math-posp y)
X     (math-add-float (math-arctan-raw (math-div-float y x))
X             (math-pi)))
X    (t
X     (math-sub-float (math-arctan-raw (math-div-float y x))
X             (math-pi))))
X)
X
X(defun math-arc-sin-cos (x)   ; [V N] [Public]
X  (if (and (Math-vectorp x)
X       (= (length x) 3))
X      (math-arctan2 (nth 2 x) (nth 1 x))
X    (math-reject-arg x "Two-element vector expected"))
X)
X(fset 'calcFunc-arcsincos (symbol-function 'math-arc-sin-cos))
X
X
X
X;;; Exponential function.
X
X(defun math-exp (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2 (math-exp-raw (math-float x))))
X    ((eq (car-safe x) 'sdev)
X     (let ((ex (math-exp (nth 1 x))))
X       (math-make-sdev ex (math-mul (nth 2 x) ex))))
X    ((eq (car-safe x) 'intv)
X     (math-make-intv (nth 1 x) (math-exp (nth 2 x)) (math-exp (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-exp x)))
X)
X(fset 'calcFunc-exp (symbol-function 'math-exp))
X
X(defun math-exp-minus-1 (x)   ; [N N] [Public]
X  (cond ((math-zerop x) '(float 0 0))
X    (calc-symbolic-mode (signal 'inexact-result nil))
X    ((Math-numberp x)
X     (math-with-extra-prec 2
X       (let ((x (math-float x)))
X         (if (and (eq (car x) 'float)
X              (math-lessp-float x '(float 1 0))
X              (math-lessp-float '(float -1 0) x))
X         (math-exp-minus-1-raw x)
X           (math-add (math-exp-raw x) -1)))))
X    ((eq (car-safe x) 'sdev)
X     (if (math-constp x)
X         (let ((ex (math-exp-minus-1 (nth 1 x))))
X           (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
X       (math-make-sdev (math-exp-minus-1 (nth 1 x))
X               (math-mul (nth 2 x) (math-exp (nth 1 x))))))
X    ((eq (car-safe x) 'intv)
X     (math-make-intv (nth 1 x)
X             (math-exp-minus-1 (nth 2 x))
X             (math-exp-minus-1 (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-expm1 x)))
X)
X(fset 'calcFunc-expm1 (symbol-function 'math-exp-minus-1))
X
X(defun math-exp10 (x)   ; [N N] [Public]
X  (math-pow '(float 1 1) x)
X)
X(fset 'calcFunc-exp10 (symbol-function 'math-exp10))
X
X(defun math-exp-raw (x)   ; [N N]
X  (cond ((math-zerop x) '(float 1 0))
X    (calc-symbolic-mode (signal 'inexact-result nil))
X    ((eq (car x) 'cplx)
X     (let ((expx (math-exp-raw (nth 1 x)))
X           (sc (math-sin-cos-raw (nth 2 x))))
X       (list 'cplx
X         (math-mul-float expx (cdr sc))
X         (math-mul-float expx (car sc)))))
X    ((eq (car x) 'polar)
X     (let ((xc (math-complex x)))
X       (list 'polar
X         (math-exp-raw (nth 1 x))
X         (nth 2 x))))
X    ((or (math-lessp-float '(float 5 -1) x)
X         (math-lessp-float x '(float -5 -1)))
X     (let* ((two-x (math-mul-float x '(float 2 0)))
X        (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
X        (hfrac (math-sub-float x (math-mul-float (math-float hint)
X                             '(float 5 -1)))))
X       (math-mul-float (math-ipow (math-sqrt-e) hint)
X               (math-add-float '(float 1 0)
X                       (math-exp-minus-1-raw hfrac)))))
X    (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))
X)
X
X(defun math-exp-minus-1-raw (x)   ; [F F]
X  (math-exp-series x 2 3 x x)
X)
X
X(defun math-exp-series (sum nfac n xpow x)
X  (math-working "exp" sum)
X  (let* ((nextx (math-mul-float xpow x))
X      (nextsum (math-add-float sum (math-div-float nextx
X                               (math-float nfac)))))
X     (if (math-nearly-equal-float sum nextsum)
X     sum
X       (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))
X)
X
X
X
X;;; Logarithms.
X
X(defun math-ln (x)   ; [N N] [Public]
X  (cond ((math-zerop x)
X     (math-reject-arg x "Logarithm of zero"))
X    ((Math-numberp x)
X     (math-with-extra-prec 2 (math-ln-raw (math-float x))))
X    ((and (eq (car-safe x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (math-posp (nth 1 x))))
X     (math-make-sdev (math-ln (nth 1 x))
X             (math-div (nth 2 x) (nth 1 x))))
X    ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X     (math-make-intv (nth 1 x) (math-ln (nth 2 x)) (math-ln (nth 3 x))))
X    ((equal x '(var e var-e))
X     1)
X    ((and (eq (car-safe x) '^)
X          (equal (nth 1 x) '(var e var-e)))
X     (nth 2 x))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-ln x)))
X)
X(fset 'calcFunc-ln (symbol-function 'math-ln))
X
X(defun math-log10 (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (let ((xf (math-float x)))
X         (if (eq (nth 1 xf) 0)
X         (math-reject-arg x "Logarithm of zero"))
X         (if (Math-integer-posp (nth 1 xf))
X         (if (eq (nth 1 xf) 1)    ; log10(1*10^n) = n
X             (math-float (nth 2 xf))
X           (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
X             (math-add-float
X              (math-div-float (math-ln-raw-2
X                       (list 'float (nth 1 xf) (- xdigs)))
X                      (math-ln-10))
X              (math-float (+ (nth 2 xf) xdigs)))))
X           (math-div (math-ln xf) (math-ln-10))))))
X    ((and (eq (car-safe x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (math-posp (nth 1 x))))
X     (math-make-sdev (math-log10 (nth 1 x))
X             (math-div (nth 2 x)
X                   (math-mul (nth 1 x) (math-ln 10)))))
X    ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X     (math-make-intv (nth 1 x)
X             (math-log10 (nth 2 x))
X             (math-log10 (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-log10 x)))
X)
X(fset 'calcFunc-log10 (symbol-function 'math-log10))
X
X(defun calcFunc-pow10 (x)
X  (math-pow '(float 1 1) x)
X)
X
X(defun math-log (x b)   ; [N N N] [Public]
X  (cond ((or (eq b 10) (equal b '(float 1 1)))
X     (math-log10 x))
X    ((math-zerop x)
X     (math-reject-arg x "Logarithm of zero"))
X    ((math-zerop b)
X     (math-reject-arg b "Logarithm of zero"))
X    ((and (Math-numberp x) (Math-numberp b))
X     (math-with-extra-prec 2
X       (math-div (math-ln-raw (math-float x))
X             (math-log-base-raw b))))
X    ((and (eq (car-safe x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (math-posp (nth 1 x)))
X          (Math-numberp b))
X     (math-make-sdev (math-log (nth 1 x) b)
X             (math-div (nth 2 x)
X                   (math-mul (nth 1 x)
X                         (math-log-base-raw b)))))
X    ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X     (math-make-intv (nth 1 x)
X             (math-log (nth 2 x) b)
X             (math-log (nth 3 x) b)))
X    (t (if (Math-numberp b)
X           (calc-record-why 'numberp x)
X         (calc-record-why 'numberp b))
X       (list 'calcFunc-log x b)))
X)
X
X(defun calcFunc-log (x &optional b)
X  (if b
X      (if (or (eq b 10) (equal b '(float 1 1)))
X      (math-normalize (list 'calcFunc-log10 x))
X    (if (equal b '(var e var-e))
X        (math-normalize (list 'calcFunc-ln x))
X      (math-log x b)))
X    (math-normalize (list 'calcFunc-ln x)))
X)
X(defun calcFunc-ilog (x &optional b)
X  (if b
X      (if (equal b '(var e var-e))
X      (math-normalize (list 'calcFunc-exp x))
X    (math-pow b x))
X    (math-normalize (list 'calcFunc-exp x)))
X)
X
X
X(defun math-log-base-raw (b)   ; [N N]
X  (if (not (and (equal (car math-log-base-cache) b)
X        (eq (nth 1 math-log-base-cache) calc-internal-prec)))
X      (setq math-log-base-cache (list b calc-internal-prec
X                      (math-ln-raw (math-float b)))))
X  (nth 2 math-log-base-cache)
X)
X(setq math-log-base-cache nil)
X
X(defun math-ln-plus-1 (x)   ; [N N] [Public]
X  (cond ((Math-equal-int x -1) (math-reject-arg x "Logarithm of zero"))
X    ((math-zerop x) '(float 0 0))
X    (calc-symbolic-mode (signal 'inexact-result nil))
X    ((Math-numberp x)
X     (math-with-extra-prec 2
X       (let ((x (math-float x)))
X         (if (and (eq (car x) 'float)
X              (math-lessp-float x '(float 5 -1))
X              (math-lessp-float '(float -5 -1) x))
X         (math-ln-plus-1-raw x)
X           (math-ln-raw (math-add-float x '(float 1 0)))))))
X    ((and (eq (car-safe x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (Math-lessp -1 (nth 1 x))))
X     (math-make-sdev (math-ln-plus-1 (nth 1 x))
X             (math-div (nth 2 x) (math-add (nth 1 x) 1))))
X    ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X     (math-make-intv (nth 1 x)
X             (math-ln-plus-1 (nth 2 x))
X             (math-ln-plus-1 (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-lnp1 x)))
X)
X(fset 'calcFunc-lnp1 (symbol-function 'math-ln-plus-1))
X
X(defun math-ln-raw (x)    ; [N N] --- must be float format!
X  (cond ((eq (car-safe x) 'cplx)
X     (list 'cplx
X           (math-mul-float (math-ln-raw
X                (math-add-float (math-sqr-float (nth 1 x))
X                        (math-sqr-float (nth 2 x))))
X                   '(float 5 -1))
X           (math-arctan2-raw (nth 2 x) (nth 1 x))))
X    ((eq (car x) 'polar)
X     (math-polar (list 'cplx
X               (math-ln-raw (nth 1 x))
X               (nth 2 x))))
X    ((Math-equal-int x 1)
X     '(float 0 0))
X    (calc-symbolic-mode (signal 'inexact-result nil))
X    ((math-posp (nth 1 x))    ; positive and real
X     (let ((xdigs (1- (math-numdigs (nth 1 x)))))
X       (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
X               (math-mul-float (math-float (+ (nth 2 x) xdigs))
X                       (math-ln-10)))))
X    ((math-zerop x)
X     (error "Logarithm of zero"))
X    ((eq calc-complex-mode 'polar)    ; negative and real
X     (math-polar
X      (list 'cplx   ; negative and real
X        (math-ln-raw (math-neg-float x))
X        (math-pi))))
X    (t (list 'cplx   ; negative and real
X         (math-ln-raw (math-neg-float x))
X         (math-pi))))
X)
X
X(defun math-ln-raw-2 (x)    ; [F F]
X  (cond ((math-lessp-float '(float 14 -1) x)
X     (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
X             (math-ln-2)))
X    (t    ; now .7 < x <= 1.4
X     (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
X                    (math-add-float x '(float 1 0))))))
X)
X
X(defun math-ln-raw-3 (x)   ; [F F]
X  (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
X          '(float 2 0))
X)
X
X;;; Compute ln((1+x)/(1-x))
X(defun math-ln-raw-series (sum n x xsqr)
X  (math-working "log" sum)
X  (let* ((nextx (math-mul-float x xsqr))
X     (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X    (if (math-nearly-equal-float sum nextsum)
X    sum
X      (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))
X)
X
X(defun math-ln-plus-1-raw (x)
X  (math-lnp1-series x 2 x (math-neg x))
X)
X
X(defun math-lnp1-series (sum n xpow x)
X  (math-working "lnp1" sum)
X  (let* ((nextx (math-mul-float xpow x))
X     (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X    (if (math-nearly-equal-float sum nextsum)
X    sum
X      (math-lnp1-series nextsum (1+ n) nextx x)))
X)
X
X(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
X  (math-ln-raw-2 '(float 1 1)))
X
X(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)
X  (math-ln-raw-3 (math-float '(frac 1 3))))
X
X
X
X;;; Hyperbolic functions.
X
X(defun math-sinh (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (let ((expx (math-exp-raw (math-float x))))
X         (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
X    ((eq (car-safe x) 'sdev)
X     (math-make-sdev (math-sinh (nth 1 x))
X             (math-mul (nth 2 x) (math-cosh (nth 1 x)))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-sinh (nth 2 x))
X             (math-sinh (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-sinh x)))
X)
X(fset 'calcFunc-sinh (symbol-function 'math-sinh))
X
X(defun math-cosh (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (let ((expx (math-exp-raw (math-float x))))
X         (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
X    ((eq (car-safe x) 'sdev)
X     (math-make-sdev (math-cosh (nth 1 x))
X             (math-mul (nth 2 x)
X                   (math-abs (math-sinh (nth 1 x))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (setq x (math-abs x))
X     (math-sort-intv (nth 1 x)
X             (math-cosh (nth 2 x))
X             (math-cosh (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-cosh x)))
X)
X(fset 'calcFunc-cosh (symbol-function 'math-cosh))
X
X(defun math-tanh (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (let* ((expx (math-exp (math-float x)))
X          (expmx (math-div 1 expx)))
X         (math-div (math-sub expx expmx)
X               (math-add expx expmx)))))
X    ((eq (car-safe x) 'sdev)
X     (math-make-sdev (math-tanh (nth 1 x))
X             (math-div (nth 2 x)
X                   (math-sqr (math-cosh (nth 1 x))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-tanh (nth 2 x))
X             (math-tanh (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-tanh x)))
X)
X(fset 'calcFunc-tanh (symbol-function 'math-tanh))
X
X(defun math-arcsinh (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
X                                 '(float 1 0)))))))
X    ((eq (car-safe x) 'sdev)
X     (math-make-sdev (math-arcsinh (nth 1 x))
X             (math-div (nth 2 x)
X                   (math-sqrt
X                    (math-add (math-sqr (nth 1 x)) 1)))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-arcsinh (nth 2 x))
X             (math-arcsinh (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-arcsinh x)))
X)
X(fset 'calcFunc-arcsinh (symbol-function 'math-arcsinh))
X
X(defun math-arccosh (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (if (or t    ; need to do this even in the real case!
X           (memq (car-safe x) '(cplx polar)))
X           (let ((xp1 (math-add 1 x)))    ; this gets the branch cuts right
X         (math-ln-raw
X          (math-add x (math-mul xp1
X                    (math-sqrt-raw (math-div (math-sub
X                                  x
X                                  '(float 1 0))
X                                 xp1))))))
X         (math-ln-raw
X          (math-add x (math-sqrt-raw (math-add (math-sqr x)
X                           '(float -1 0))))))))
X    ((and (eq (car-safe x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (not (Math-lessp (nth 1 x) 1))))
X     (math-make-sdev (math-arccosh (nth 1 x))
X             (math-div (nth 2 x)
X                   (math-sqrt
X                    (math-add (math-sqr (nth 1 x)) -1)))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-arccosh (nth 2 x))
X             (math-arccosh (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-arccosh x)))
X)
X(fset 'calcFunc-arccosh (symbol-function 'math-arccosh))
X
X(defun math-arctanh (x)   ; [N N] [Public]
X  (cond ((Math-numberp x)
X     (math-with-extra-prec 2
X       (if (memq (car-safe x) '(cplx polar))
X           (math-ln-raw
X        (math-mul (math-add 1 x)
X              (math-sqrt-raw
X               (math-div '(float 1 0) (math-sub 1 (math-sqr x))))))
X         (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
X                          (math-sub 1 x)))
X               '(float 5 -1)))))
X    ((and (eq (car-safe x) 'sdev)
X          (or (not (math-constp (nth 1 x)))
X          (Math-lessp (math-abs (nth 1 x)) 1)))
X     (math-make-sdev (math-arctanh (nth 1 x))
X             (math-div (nth 2 x)
X                   (math-sub 1 (math-sqr (nth 1 x))))))
X    ((and (eq (car x) 'intv) (math-constp x))
X     (math-sort-intv (nth 1 x)
X             (math-arctanh (nth 2 x))
X             (math-arctanh (nth 3 x))))
X    (t (calc-record-why 'numberp x)
X       (list 'calcFunc-arctanh x)))
X)
X(fset 'calcFunc-arctanh (symbol-function 'math-arctanh))
X
X
X;;; Convert A from HMS or degrees to radians.
X(defun math-deg-to-rad (a)   ; [R R] [Public]
X  (cond ((or (Math-numberp a)
X         (eq (car a) 'intv))
X     (math-mul a (math-pi-over-180)))
X    ((eq (car a) 'hms)
X     (math-from-hms a 'rad))
X    ((eq (car a) 'sdev)
X     (math-make-sdev (math-deg-to-rad (nth 1 a))
X             (math-deg-to-rad (nth 2 a))))
X    (t (list 'calcFunc-rad a)))
X)
X(fset 'calcFunc-rad (symbol-function 'math-deg-to-rad))
X
X;;; Convert A from HMS or radians to degrees.
X(defun math-rad-to-deg (a)   ; [R R] [Public]
X  (cond ((or (Math-numberp a)
X         (eq (car a) 'intv))
X     (math-div a (math-pi-over-180)))
X    ((eq (car a) 'hms)
X     (math-from-hms a 'deg))
X    ((eq (car a) 'sdev)
X     (math-make-sdev (math-rad-to-deg (nth 1 a))
X             (math-rad-to-deg (nth 2 a))))
X    (t (list 'calcFunc-deg a)))
X)
X(fset 'calcFunc-deg (symbol-function 'math-rad-to-deg))
X
X
X
X
X;;;; Number theory.
X
X(defun calcFunc-idiv (a b)   ; [I I I] [Public]
X  (cond ((and (Math-natnump a) (Math-natnump b) (not (eq b 0)))
X     (math-quotient a b))
X    ((Math-realp a)
X     (if (Math-realp b)
X         (let ((calc-prefer-frac t))
X           (math-floor (math-div a b)))
X       (math-reject-arg b 'realp)))
X    ((eq (car-safe a) 'hms)
X     (if (eq (car-safe b) 'hms)
X         (let ((calc-prefer-frac t))
X           (math-floor (math-div a b)))
X       (math-reject-arg b 'hmsp)))
X    ((and (or (eq (car-safe a) 'intv) (Math-realp a))
X          (or (eq (car-safe b) 'intv) (Math-realp b)))
X     (math-floor (math-div a b)))
X    (t (math-reject-arg a 'anglep)))
X)
X
X(defun calcFunc-fdiv (a b)   ; [R I I] [Public]
X  (if (Math-num-integerp a)
X      (if (Math-num-integerp b)
X      (if (Math-zerop b)
X          (math-reject-arg a "Division by zero")
X        (math-make-frac (math-trunc a) (math-trunc b)))
X    (math-reject-arg b 'integerp))
X    (math-reject-arg a 'integerp))
X)
X
X(defun math-lcm (a b)
X  (let ((g (math-gcd a b)))
X    (if (Math-numberp g)
X    (math-div (math-mul a b) g)
X      (list 'calcFunc-lcm a b)))
X)
X(fset 'calcFunc-lcm (symbol-function 'math-lcm))
X
X(defun math-extended-gcd (a b)   ; Knuth section 4.5.2
X  (cond
X   ((not (Math-integerp a))
X    (if (Math-messy-integerp a)
X    (math-extended-gcd (math-trunc a) b)
X      (calc-record-why 'integerp a)
X      (list 'calcFunc-egcd a b)))
X   ((not (Math-integerp b))
X    (if (Math-messy-integerp b)
X    (math-extended-gcd a (math-trunc b))
X      (calc-record-why 'integerp b)
X      (list 'calcFunc-egcd a b)))
X   (t
X    (let ((u1 1) (u2 0) (u3 a)
X      (v1 0) (v2 1) (v3 b)
X      t1 t2 q)
X      (while (not (eq v3 0))
X    (setq q (math-idivmod u3 v3)
X          t1 (math-sub u1 (math-mul v1 (car q)))
X          t2 (math-sub u2 (math-mul v2 (car q)))
X          u1 v1  u2 v2  u3 v3
X          v1 t1  v2 t2  v3 (cdr q)))
X      (list 'vec u3 u1 u2))))
X)
X(fset 'calcFunc-egcd (symbol-function 'math-extended-gcd))
X
X
X;;; Factorial and related functions.
X
X(defun math-factorial (n)   ; [I I] [F F] [Public]
X  (let (temp)
X    (cond ((Math-integer-negp n) (list 'calcFunc-fact n))
X      ((Math-zerop n) 1)
X      ((integerp n) (math-factorial-iter (1- n) 2 1))
X      ((and (math-messy-integerp n)
X        (Math-lessp (setq temp (math-trunc n)) 100))
X       (if (natnump temp)
X           (math-with-extra-prec 1
X         (math-factorial-iter (1- temp) 2 '(float 1 0)))
X         (list 'calcFunc-fact max)))
X      ((math-realp n)
X       (math-with-extra-prec 3
X         (math-gammap1-raw (math-float n))))
X      (t (calc-record-why 'realp n)
X         (list 'calcFunc-fact n))))
X)
X(fset 'calcFunc-fact (symbol-function 'math-factorial))
X
X(defun math-factorial-iter (count n f)
X  (if (= (% n 5) 1)
X      (math-working (format "factorial(%d)" (1- n)) f))
X  (if (> count 0)
X      (math-factorial-iter (1- count) (1+ n) (math-mul n f))
X    f)
X)
X
X(math-defcache math-sqrt-two-pi nil
X  (math-sqrt (math-two-pi)))
X
X(defun math-gammap1-raw (x)   ; compute gamma(1 + x)
X  (cond ((math-lessp-float x '(float 1 1))
X     (if (math-lessp-float x '(float -10 0))
X         (setq x (math-neg-float
X              (math-div-float
X               (math-pi)
X               (math-mul-float (math-gammap1-raw
X                    (math-add-float (math-neg-float x)
X                            '(float -1 0)))
X                       (math-sin-raw
X                    (math-mul (math-pi) x)))))))
X     (let ((xplus1 (math-add-float x '(float 1 0))))
X       (math-div-float (math-gammap1-raw xplus1) xplus1)))
X    (t   ; x now >= 10.0
X     (let ((xinv (math-div 1 x))
X           (lnx (math-ln-raw x)))
X       (math-mul (math-sqrt-two-pi)
X             (math-exp-raw
X              (math-gamma-series
X               (math-sub (math-mul (math-add x '(float 5 -1))
X                       lnx)
X                 x)
X               xinv
X               (math-sqr xinv)
X               2))))))
X)
X
X(defun calcFunc-gamma (x)
X  (calcFunc-fact (math-add x -1))
X)
X
X(defun math-gamma-series (sum x xinvsqr n)
X  (math-working "gamma" sum)
X  (let* ((bn (math-bernoulli-number n))   ; this will always be a "frac" form.
X     (next (math-add-float
X        sum
X        (math-mul-float (math-div-float (math-float (nth 1 bn))
X                        (math-float (* (nth 2 bn)
X                                   (* n (1- n)))))
X                x))))
X    (if (math-nearly-equal-float sum next)
X    next
X      (if (= n 24)
X      (progn
X        (calc-record-why "Gamma computation stopped early, not all digits may be valid")
X        next)
X    (math-gamma-series next (math-mul-float x xinvsqr) xinvsqr (+ n 2)))))
X)
X
X(defun math-bernoulli-number (n)
X  (if (= n 1)
X      '(frac -1 2)
X    (if (= (% n 2) 1)
X    0
X      (aref '[ 1 (frac 1 6) (frac -1 30) (frac 1 42) (frac -1 30)
X           (frac 5 66) (frac -691 2730) (frac 7 6) (frac -3617 510)
X           (frac 43867 798) (frac -174611 330) (frac 854513 138)
X           (frac (bigneg 91 364 236) 2730) ]
X        (/ n 2))))
X)
X;;; To come up with more, we could use this rule:
X;;;   Bn = n! bn
X;;;   bn = - sum_k=0^n-1 bk / (n-k+1)!
X
X(defun math-double-factorial (n)   ; [I I] [F F] [Public]
X  (cond ((Math-integer-negp n) (list 'calcFunc-dfact n))
X    ((Math-zerop n) 1)
X    ((integerp n) (math-double-factorial-iter n (+ 2 (% n 2)) 1))
X    ((math-messy-integerp n)
X     (let ((temp (math-trunc n)))
X       (if (natnump temp)
X           (math-with-extra-prec 1
X         (math-double-factorial-iter temp (+ 2 (% temp 2))
X                         '(float 1 0)))
X         (list 'calcFunc-dfact max))))
X    (t (calc-record-why 'natnump n)
X       (list 'calcFunc-dfact n)))
X)
X(fset 'calcFunc-dfact (symbol-function 'math-double-factorial))
X
X(defun math-double-factorial-iter (max n f)
X  (if (< (% n 10) 2)
X      (math-working (format "dfact(%d)" (- n 2)) f))
X  (if (<= n max)
X      (math-double-factorial-iter max (+ n 2) (math-mul n f))
X    f)
X)
X
X(defun math-permutations (n m)   ; [I I I] [F F F] [Public]
X  (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
X     (math-factorial-iter n (1+ (- n m)) 1))
X    ((or (not (math-num-integerp n))
X         (not (math-num-integerp m)))
X     (or (math-realp n) (math-reject-arg n 'realp))
X     (or (math-realp m) (math-reject-arg m 'realp))
X     (and (math-num-integerp n) (math-negp n) (math-reject-arg n 'range))
X     (and (math-num-integerp m) (math-negp m) (math-reject-arg m 'range))
X     (math-div (math-factorial n) (math-factorial m)))
X    (t
X     (let ((tn (math-trunc n))
X           (tm (math-trunc m)))
X       (or (integerp tn) (math-reject-arg tn 'fixnump))
X       (or (integerp tm) (math-reject-arg tm 'fixnump))
X       (or (and (<= tm tn) (>= tm 0)) (math-reject-arg tm 'range))
X       (math-with-extra-prec 1
X         (math-factorial-iter tm (1+ (- tn tm)) '(float 1 0))))))
X)
X(fset 'calcFunc-perm (symbol-function 'math-permutations))
X
X(defun math-choose (n m)   ; [I I I] [F F F] [Public]
X  (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
X     (math-choose-iter m n 1 1))
X    ((not (math-realp n))
X     (math-reject-arg n 'realp))
X    ((not (math-realp m))
X     (math-reject-arg m 'realp))
X    ((not (math-num-integerp m))
X     (if (and (math-num-integerp n) (math-negp n))
X         (list 'calcFunc-choose n m)
X       (math-div (math-factorial (math-float n))
X             (math-mul (math-factorial m)
X                   (math-factorial (math-sub n m))))))
X    ((math-negp m) 0)
X    ((math-negp n)
X     (let ((val (math-choose (math-add (math-add n m) -1) m)))
X       (if (math-evenp (math-trunc m))
X           val
X         (math-neg val))))
X    ((and (math-num-integerp n)
X          (Math-lessp n m))
X     0)
X    (t
X     (let ((tm (math-trunc m)))
X       (or (integerp tm) (math-reject-arg tm 'fixnump))
X       (if (> tm 100)
X           (math-div (math-factorial (math-float n))
X             (math-mul (math-factorial (math-float m))
X                   (math-factorial (math-float
X                            (math-sub n m)))))
X         (math-with-extra-prec 1
X           (math-choose-float-iter tm n 1 '(float 1 0)))))))
X)
X(fset 'calcFunc-choose (symbol-function 'math-choose))
X
X(defun math-choose-iter (m n i c)
X  (if (= (% i 5) 1)
X      (math-working (format "choose(%d)" (1- i)) c))
X  (if (<= i m)
X      (math-choose-iter m (1- n) (1+ i)
X            (math-quotient (math-mul c n) i))
X    c)
X)
X
X(defun math-choose-float-iter (count n i c)
X  (if (= (% i 5) 1)
X      (math-working (format "choose(%d)" (1- i)) c))
X  (if (> count 0)
X      (math-choose-float-iter (1- count) (math-sub n 1) (1+ i)
X                  (math-div (math-mul c n) i))
X    c)
X)
X
X
X;;; Produce a random digit in the range 0..999.
X;;; Avoid various pitfalls that may lurk in the built-in (random) function!
X(defun math-random-digit ()
X  (prog1
X      (% (lsh (random math-first-random-flag) -4) 1000)
X    (setq math-first-random-flag nil))
X)
X(setq math-first-random-flag t)
X
X;;; Produce an N-digit random integer.
X(defun math-random-digits (n)
X  (cond ((<= n 6)
X     (math-scale-right (+ (* (math-random-digit) 1000) (math-random-digit))
X               (- 6 n)))
X    (t (let* ((slop (% (- 900003 n) 3))
X          (i (/ (+ n slop) 3))
X          (digs nil))
X         (while (> i 0)
X           (setq digs (cons (math-random-digit) digs)
X             i (1- i)))
X         (math-normalize (math-scale-right (cons 'bigpos digs)
X                           slop)))))
X)
X
X;;; Produce a uniformly-distributed random float 0 <= N < 1.
X(defun math-random-float ()
X  (math-make-float (math-random-digits calc-internal-prec)
X           (- calc-internal-prec))
X)
X
X;;; Produce a Gaussian-distributed random float with mean=0, sigma=1.
X(defun math-gaussian-float ()
X  (math-with-extra-prec 2
X    (if (and math-gaussian-cache
X         (= (car math-gaussian-cache) calc-internal-prec))
X    (prog1
X        (cdr math-gaussian-cache)
X      (setq math-gaussian-cache nil))
X      (let* ((v1 (math-add (math-mul (math-random-float) 2) -1))
X         (v2 (math-add (math-mul (math-random-float) 2) -1))
X         (r (math-add (math-sqr v1) (math-sqr v2))))
X    (while (or (not (Math-lessp r 1)) (math-zerop r))
X      (setq v1 (math-add (math-mul (math-random-float) 2) -1)
X        v2 (math-add (math-mul (math-random-float) 2) -1)
X        r (math-add (math-sqr v1) (math-sqr v2))))
X    (let ((fac (math-sqrt (math-mul (math-div (math-ln r) r) -2))))
X      (setq math-gaussian-cache (cons calc-internal-prec
X                      (math-mul v1 fac)))
X      (math-mul v2 fac)))))
X)
X(setq math-gaussian-cache nil)
X
X;;; Produce a random integer or real 0 <= N < MAX.
X(defun math-random (max)
X  (cond ((Math-zerop max)
X     (math-gaussian-float))
X    ((Math-integerp max)
X     (let* ((digs (math-numdigs max))
X        (r (math-random-digits (+ digs 3))))
X       (math-mod r max)))
X    ((Math-realp max)
X     (math-mul (math-random-float) max))
X    ((and (eq (car max) 'intv) (math-constp max)
X          (Math-lessp (nth 2 max) (nth 3 max)))
X     (if (math-floatp max)
X         (let ((val (math-add (math-mul (math-random-float)
X                        (math-sub (nth 3 max) (nth 2 max)))
X                  (nth 2 max))))
X           (if (or (and (memq (nth 1 max) '(0 1))      ; almost not worth
X                (math-equal val (nth 2 max)))  ;   checking!
X               (and (memq (nth 1 max) '(0 2))
X                (math-equal val (nth 3 max))))
X           (math-random max)
X         val))
X       (let ((lo (if (memq (nth 1 max) '(0 1))
X             (math-add (nth 2 max) 1) (nth 2 max)))
X         (hi (if (memq (nth 1 max) '(1 3))
X             (math-add (nth 3 max) 1) (nth 3 max))))
X         (if (Math-lessp lo hi)
X         (math-add (math-random (math-sub hi lo)) lo)
X           (math-reject-arg max "Empty interval")))))
X    (t (math-reject-arg max 'realp)))
X)
X(fset 'calcFunc-random (symbol-function 'math-random))
X
X
X;;; Check if the integer N is prime.  [X I]
X;;; Return (nil) if non-prime,
X;;;        (nil N) if non-prime with known factor N,
X;;;        (nil unknown) if non-prime with no known factors,
X;;;        (t) if prime,
X;;;        (maybe N P) if probably prime (after N iters with probability P%)
X(defun math-prime-test (n iters)
X  (if (and (Math-vectorp n) (cdr n))
X      (setq n (nth (1- (length n)) n)))
X  (if (Math-messy-integerp n)
X      (setq n (math-trunc n)))
X  (let ((res))
X    (while (> iters 0)
X      (setq res
X        (cond ((and (integerp n) (<= n 5003))
X           (list (= (math-next-small-prime n) n)))
X          ((not (Math-integerp n))
X           (error "Argument must be an integer"))
X          ((Math-integer-negp n)
X           '(nil))
X          ((Math-natnum-lessp n '(bigpos 0 0 8))
X           (setq n (math-fixnum n))
X           (let ((i -1) v)
X             (while (and (> (% n (setq v (aref math-primes-table
X                               (setq i (1+ i)))))
X                    0)
X                 (< (* v v) n)))
X             (if (= (% n v) 0)
X             (list nil v)
X               '(t))))
X          ((not (equal n (car math-prime-test-cache)))
X           (cond ((= (% (nth 1 n) 2) 0) '(nil 2))
X             ((= (% (nth 1 n) 5) 0) '(nil 5))
X             (t (let ((dig (cdr n)) (sum 0))
X                  (while dig
X                (if (cdr dig)
X                    (setq sum (% (+ (+ sum (car dig))
X                            (* (nth 1 dig) 1000))
X                         111111)
X                      dig (cdr (cdr dig)))
X                  (setq sum (% (+ sum (car dig)) 111111)
X                    dig nil)))
X                  (cond ((= (% sum 3) 0) '(nil 3))
X                    ((= (% sum 7) 0) '(nil 7))
X                    ((= (% sum 11) 0) '(nil 11))
X                    ((= (% sum 13) 0) '(nil 13))
X                    ((= (% sum 37) 0) '(nil 37))
X                    (t
X                     (setq math-prime-test-cache-k 1
X                       math-prime-test-cache-q
X                       (math-div2 n)
X                       math-prime-test-cache-nm1
X                       (math-add n -1))
X                     (while (math-evenp
X                         math-prime-test-cache-q)
X                       (setq math-prime-test-cache-k
X                         (1+ math-prime-test-cache-k)
X                         math-prime-test-cache-q
X                         (math-div2
X                          math-prime-test-cache-q)))
X                     (setq iters (1+ iters))
X                     (list 'maybe
X                       0
X                       (math-sub
X                        100
X                        (math-div
X                         '(float 232 0)
X                         (math-numdigs n))))))))))
X          ((not (eq (car (nth 1 math-prime-test-cache)) 'maybe))
X           (nth 1 math-prime-test-cache))
X          (t   ; Fermat step
X           (let* ((x (math-add (math-random (math-add n -2)) 2))
X              (y (math-pow-mod x math-prime-test-cache-q n))
X              (j 0))
X             (while (and (not (eq y 1))
X                 (not (equal y math-prime-test-cache-nm1))
X                 (< (setq j (1+ j)) math-prime-test-cache-k))
X               (setq y (math-mod (math-mul y y) n)))
X             (if (or (equal y math-prime-test-cache-nm1)
X                 (and (eq y 1) (eq j 0)))
X             (list 'maybe
X                   (1+ (nth 1 (nth 1 math-prime-test-cache)))
X                   (math-mul (nth 2 (nth 1 math-prime-test-cache))
X                     '(float 25 -2)))
X               '(nil unknown))))))
X      (setq math-prime-test-cache (list n res)
X        iters (if (eq (car res) 'maybe)
X              (1- iters)
X            0)))
X    res)
X)
X(defvar math-prime-test-cache '(-1))
X
X;;; Theory: summing base-10^6 digits modulo 111111 is "casting out 999999s".
X;;; Initial probability that N is prime is 1/ln(N) = log10(e)/log10(N).
X;;; After culling [2,3,5,7,11,13,37], probability of primality is 5.36 x more.
X;;; Initial reported probability of non-primality is thus 100% - this.
X;;; Each Fermat step multiplies this probability by 25%.
X;;; The Fermat step is algorithm P from Knuth section 4.5.4.
X
X
X(defun math-prime-factors (n)
X  (setq math-prime-factors-finished t)
X  (if (Math-messy-integerp n)
X      (setq n (math-trunc n)))
X  (if (and (Math-natnump n) (Math-natnum-lessp 2 n))
X      (let (factors res p (i 0))
X    (while (and (not (eq n 1))
X            (< i (length math-primes-table)))
X      (setq p (aref math-primes-table i))
X      (while (eq (cdr (setq res (cond ((eq n p) (cons 1 0))
X                      ((eq n 1) (cons 0 1))
X                      ((consp n) (math-idivmod n p))
X                      (t (cons (/ n p) (% n p))))))
X             0)
X        (math-working "factor" p)
X        (setq factors (nconc factors (list p))
X          n (car res)))
X      (or (eq n 1)
X          (Math-natnum-lessp p (car res))
X          (setq factors (nconc factors (list n))
X            n 1))
X      (setq i (1+ i)))
X    (or (setq math-prime-factors-finished (eq n 1))
X        (setq factors (nconc factors (list n))))
X    (cons 'vec factors))
X    (calc-record-why 'integerp n)
X    (list 'calcFunc-prfac n))
X)
X(fset 'calcFunc-prfac (symbol-function 'math-prime-factors))
X
X(defun math-totient (n)
X  (if (Math-messy-integerp n)
X      (setq n (math-trunc n)))
X  (if (Math-natnump n)
X      (if (Math-natnum-lessp n 2)
X      (if (Math-negp n)
X          (math-totient (math-abs n))
X        n)
X    (let ((factors (cdr (math-prime-factors n)))
X          p)
X      (if math-prime-factors-finished
X          (progn
X        (while factors
X          (setq p (car factors)
X            n (math-mul (math-div n p) (math-add p -1)))
X          (while (equal p (car factors))
X            (setq factors (cdr factors))))
X        n)
X        (calc-record-why "Number too big to factor" n)
X        (list 'calcFunc-totient n))))
X    (calc-record-why 'natnump n)
X    (list 'calcFunc-totient n))
X)
X(fset 'calcFunc-totient (symbol-function 'math-totient))
X
X(defun math-moebius (n)
X  (if (Math-messy-integerp n)
X      (setq n (math-trunc n)))
X  (if (and (Math-natnump n) (not (eq n 0)))
X      (if (Math-natnum-lessp n 2)
X      (if (Math-negp n)
X          (math-moebius (math-abs n))
X        1)
X    (let ((factors (cdr (math-prime-factors n)))
X          (mu 1))
X      (if math-prime-factors-finished
X          (progn
X        (while factors
X          (setq mu (if (equal (car factors) (nth 1 factors))
X                   0 (math-neg mu))
X            factors (cdr factors)))
X        mu)
X        (calc-record-why "Number too big to factor" n)
X        (list 'calcFunc-moebius n))))
X    (calc-record-why 'natnump n)
X    (list 'calcFunc-moebius n))
X)
X(fset 'calcFunc-moebius (symbol-function 'math-moebius))
X
X
X(defun math-next-prime (n iters)
X  (if (Math-integerp n)
X      (if (Math-integer-negp n)
X      2
X    (if (and (integerp n) (< n 5003))
X        (math-next-small-prime (1+ n))
X      (if (math-evenp n)
X          (setq n (math-add n -1)))
X      (let (res)
X        (while (not (car (setq res (math-prime-test
X                    (setq n (math-add n 2)) iters)))))
X        (if (and calc-verbose-nextprime
X             (eq (car res) 'maybe))
X        (calc-report-prime-test res)))
X      n))
X    (if (Math-realp n)
X    (math-next-prime (math-trunc n) iters)
X      (math-reject-arg n 'integerp)))
X)
X(fset 'calcFunc-nextprime (symbol-function 'math-next-prime))
X(setq calc-verbose-nextprime nil)
X
X(defun math-prev-prime (n iters)
X  (if (Math-integerp n)
X      (if (Math-lessp n 4)
X      2
X    (if (math-evenp n)
X        (setq n (math-add n 1)))
X    (let (res)
X      (while (not (car (setq res (math-prime-test
X                      (setq n (math-add n -2)) iters)))))
X      (if (and calc-verbose-nextprime
X           (eq (car res) 'maybe))
X          (calc-report-prime-test res)))
X    n)
X    (if (Math-realp n)
X    (math-prev-prime (math-ceiling n) iters)
X      (math-reject-arg n 'integerp)))
X)
X(fset 'calcFunc-prevprime (symbol-function 'math-prev-prime))
X
X(defun math-next-small-prime (n)
X  (if (and (integerp n) (> n 2))
X      (let ((lo -1)
X        (hi (length math-primes-table))
X        mid)
X    (while (> (- hi lo) 1)
X      (if (> n (aref math-primes-table
X             (setq mid (ash (+ lo hi) -1))))
X          (setq lo mid)
X        (setq hi mid)))
X    (aref math-primes-table hi))
X    2)
X)
X
X(defconst math-primes-table
X  [2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89
X     97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181
X     191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277
X     281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383
X     389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487
X     491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601
X     607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709
X     719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827
X     829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947
X     953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049
X     1051 1061 1063 1069 1087 1091 1093 1097 1103 1109 1117 1123 1129 1151
SHAR_EOF
echo "End of part 8"
echo "File calc-ext.el is continued in part 9"
echo "9" > s2_seq_.tmp
exit 0